Molecular dynamics of shock fronts and their transitions in Lennard-Jonesium and Tin 
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We develop a Continuous Hugoniot Method for the efficient simulation of shock wave fronts with 
molecular dynamics. This approach achieves a significantly improved efficiency in the generation of 
a dense sampling of steady-state shock front states, and allows for the study of shocks as a function 
of a continuous shock strength parameter, v p . This is, to our knowledge, the first attempt to map 
out the Hugoniot in a continuous fashion. 

We first apply this method to shocks in single-crystal Lennard-Jonesium along the (f00) direc- 
tion. Excellent agreement is found with both the published Lennard- Jones Hugoniot and results of 
conventional simulation methods. 

We next present a continuous numerical Hugoniot for shocks in tin which agrees to within 6% 
with experimental data. We study the strong shock to elastic-plastic shock transition in tin and 
find that it is a continuous transition consistent with a transcritical bifurcation. 
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I. INTRODUCTION 

Laser-induced shock waves and femtosecond time- 
resolved optical diagnostics are providing a new win- 
dow into shock research. Recent shock physics exper- 
iments have made atomic response measurements at 
length scales and time scales well matched with computa- 
tional molecular dynamics. There are new opportunities 
for direct collaboration between experiment and simula- 
tion in the study of dynamics at the shock front. 

This work presents a method which allows effective 
collaboration with experiment in characterizing the dy- 
namics and melt properties at the shock front in real 
materials. Recent experiments have indicated that the 
dynamic response at the shock front in solids can depend 
strongly not only on shock strength but also on failure 
or transition timescales . On-going experiments at the 
High-Intensity Laser Science Group at the University of 
Texas at Austin are investigating the shock-strength de- 
pendence of melt time scales in tin and aluminum using 
table-top and larger laser systems in the terawatt range. 

There are two major difficulties in applying conven- 
tional simulation methods to the study of dynamics near 
a shock front: (1) To produce the environment at the 
front, one must simulate a large and ever-growing system, 
of which the front constitutes only a very small fraction; 
and (2) The conditions within a steady-state shock take 
long times to arise, and each computationally-expensive 
shock run results in only a single data point. 

Our goal in this paper is to address each of these de- 
ficiencies with an efficient approach. Our method, in- 
troduced in Section II, focuses computational resources 
only on the shock front, while simultaneously producing 
a continuum of shock strength final states in a single run. 

The constrained dynamics methods of the Hugoniostat 
and others offer solutions to the first problem, but 



offer no information about the non-equilibrium dynamics 
at the shock front. The approach of Zhakhovskii et al. 
[f| succeeds in addressing the first point at the shock 
front, but does not address the second. We generalize 
and expand on these methods here. First, we concentrate 
our efforts on the neighborhood of the shock interface, 
thereby increasing computational efficiency, and second, 
we map system response to a continuum of shock strength 
final states by continuously varying the parameter within 
the simulation. This combination of techniques we call 
the Continuous Hugoniot Method. 

The shock Hugoniot is a relationship between any two 
variables of the final material state behind a shock. It is 
derived from the conservation of mass, momentum, and 
energy across a shock front and the equation of state 
(EOS) for the material. The conservation equations are 
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where U, u, p — y, P and E, are the shock velocity, par- 
ticle velocity, density, pressure tensor and internal energy, 
respectively. The subscript indicates the initial mate- 
rial state ahead of the shock front. The state variables 
refer to values far behind the shock front. P xx refers 
to the component of the pressure tensor far behind the 
shock front perpendicular to the direction of the front. 
The EOS is material specific. The Hugoniot, a curve de- 
rived from these four equations and five unknowns, is a 
locus of final states. Each point on the Hugoniot curve 
represents the final state of an individual shock experi- 
ment or simulation. 
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II. CONTINUOUS HUGONIOT METHOD 

Introductions to shock physics Q state clearly that 
the Hugoniot is not a thermodynamic path. However, our 
simulation method makes the Hugoniot the thermody- 
namic path of our simulation system. This is not thought 
to be possible experimentally. Figure shows (A) the ex- 
perimental loading paths (Rayleigh lines) to each state of 
the Hugoniot and (B) the loading path which we will use 
to reach the same state points. The following was briefly 
summarized in [? ]. 




Volume Volume 

FIG. 1: Hugoniot as a destination versus path - (A) 

The Hugoniot (dotted line) is illustrated as a collection of 
final shock states. (B) The Continuous Hugoniot Method 
follows the path of the Hugoniot continuously from the initial 
reference state. The solid lines in both plots depict the path 
of the system through P-V space. 



Conventional Simulation Reference Runs 

The Continuous Hugoniot Method results in a contin- 
uum of shock states which range from some initial shock 
strength to some final shock strength. Two conventional 
molecular dynamics shock wave computations bookend 
this continuum. Conventional simulations are runs which 
arc driven by a warm impactor (momentum mirror) and 
are allowed to evolve for long times until they converge 
to a steady state. These reference systems generally grow 
very large and are therefore computationally expensive. 
Runs of this type are, however, necessary to seed our 
method and bound its error. 

The method begins with a conventional run which 
leaves the system in the initial Hugoniot reference state. 
This run is illustrated in Fig. [21 (B) by the straight line 
chord (Rayleigh Line) connecting the initial condition to 
the Hugoniot. This reference run provides a initial Hugo- 
niot state. We determine the smallest system necessary 
to model the shock front from this state. The reduced 
system size is found by measuring the distance behind 
the shock front at which the material has thermalized. 



Thermalization is determined from velocity distribution 
analysis at increasing distances from the shock front. For 
our purposes, we consider the system thermalized when 
the distributions in the transverse and propogation di- 
rections are Maxwell-Boltzmann and give the same tem- 
peratures. 

System Reduction 

We truncate the reference system to a fixed width, leav- 
ing pristine material ahead and slab of compressed mate- 
rial behind. The reduced system width is determined by 
the reference state thermalization length we have mea- 
sured. This reduced system is then evolved with ma- 
terial added and purged based on criterion established 
to keep the shock front in the center of the fixed-width 
sample. At the rear purge point, we preserve the velocity 
distributions with a Langevin thermostat, and a momen- 
tum mirror prevents particle loss. The piston velocity 
and mean velocity of the thermostat are equal, and this 
value, Up, serves as the external control parameter for 
shock strength. 

A buffer zone of undisturbed crystal must always ex- 
ist ahead of the shock front, but this buffer need not be 
large. The shock travels faster than any wave in the un- 
compressed material, so we do not need to be concerned 
with preheating or dispersion. The buffer zone we main- 
tain is 5 to 10 lattice planes wide. The initial tempera- 
ture is set by another Langevin thermostat, and particles 
within a unit cell of the system's forward edge are frozen 
in place to prevent unphysical surface reconstructions. 

As the system evolves, the shock front consumes ma- 
terial and advances into the buffer zone. The shock 
front location is continually updated at a set interval of 
timesteps. The front is determined by analysis of the 
gradient in total energy per particle along the propoga- 
tion direction. The front location is used to trigger an 
advance event when the size of the buffer zone (measured 
from the system's forward edge to the shock front) falls 
below the user defined minimum. 

At the advance event, a block of crystal is generated 
and appended to the forward system boundary. The gen- 
erated crystal must have thickness equal to an integer 
multiple of the crystal unit cell and have zero initial tem- 
perature to prevent a mismatch with the frozen forward 
edge of the system. 

The requirements for a smooth and valid purge event 
will be discussed at end of this section. For now, we as- 
sume the requirements are satisfied and discuss the steps 
of the purge process. At each advance event, the total 
system size is calculated. A purge is initiated if the sys- 
tem size exceeds a user-defined maximum. The trunca- 
tion of the system behind the front can be very delicate 
because the system is driven from behind. Care must 
be taken to preserve the thermodynamic properties at 
the back edge of the system and not to introduce any 
unphysical disturbance that might propagate. 

Material is purged in a block equal in volume to the 
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FIG. 2: Co-moving frame centered on the shock front - From an initial reference state (A) a reduced system of particles 
are extracted. A thermostated piston drives the system from the left and particles on the right are frozen. (B) The system 
evolves as the piston drives forward. As the shock front advances, fresh material is added to the right. (C) When the material 
at the left has thermalized ahead of the piston, a block of material is purged and the cycle repeats. 



block appended during an advance. Because the purged 
material is compressed, however, there are more particles 
purged than appended in a single event. For this rea- 
son, there are fewer purge events than advance events. 
A purge begins with a temperature measurement at the 
purge point. The back boundary is then shifted forward 
and particles beyond the new boundary are dropped from 
the system arrays. From this point, the system is evolved 
normally. 

Our goals in this system reduction method are to sim- 
ulate without approximation the dynamics at the front 
of a shock wave without devoting computing resources to 
simulate the entire system behind the shock. To do this, 
we require (1) that the final shock be in a steady state, 
and (2) that a thermalized equilibrium thermodynamic 
state evolve quickly (i.e. within a short distance) behind 
the shock front. In some cases, this thermodynamic state 
is the equilibrium final state given by the Hugoniot rela- 
tion. In others, it is an intermediate state. Our method 
requires only that we can replicate the thermodynamic 
properties at the purge with a Maxwell-Boltzmann ve- 
locity distribution. 

It was observed early in this work that the piston 
boundary conditions play a very important role in the 
ability to drive the system smoothly. It must be re- 
membered here that the Shockwave is supersonic in the 
uncompressed material, but subsonic in the compressed 
material behind the front. This means that disturbances 
from the piston travel faster than the shock front and 
will eventually interact with it. In traditional shock sim- 



ulations, one does not need to worry about the piston 
interaction because as the piston grows further and fur- 
ther isolated from the shock front the two become de- 
coupled and interact only in an average thermodynamic 
way. In these large systems, any piston-generated spatial 
correlations in the system decay away in a trivially small 
percentage of the whole system. 

In our reduced systems, we need to be much more 
careful with the piston interaction. The piston can- 
not be allowed to produce spatial correlation, and the 
purge events must be smooth to ensure that no periodic 
pulses are introduced to disturb the front dynamics. Fig- 
ure [S] demonstrates the problems generated by driving a 
tin melt with a conventional momentum mirror alone. 
Driving the same system with a momentum mirror and 
strong stochastic forcing region is also shown. The spa- 
tial correlations are clearly visible without the stochastic 
forcing (top) and disappear when it is introduced (bot- 
tom). Weak Langevin thermostats showed results simi- 
lar to those of the momentum mirror alone. We found 
that damping coefficients very nearly equal to 1 /timestep 
were necessary to prevent the effect. For these studies 
the Langevin thermostat was implemented with particle 
damping and random kicking. The amplitude of the kick- 
ing was set by the system parameters, using the modified 
equation of motion, 



F = ma — mbv + £r (4) 
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FIG. 3: Piston-induced spatial correlations - Two par- 
ticle slices illustrate the potential problem of spatial correla- 
tions developing near the piston boundary, (top) The mo- 
mentum mirror alone shows a clear particle concentration at 
the piston surface with a first-neighbor gap clearly develop- 
ing, (bottom) Adding a region of stochastic forcing ahead of 
the momentum mirror removes the problem. 



where 
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(5) 



is the amplitude of the stochastic force and r is a random- 
ized dimensionless vector with components ranging from 
-1 to 1. Here m is the mass, b is the damping coefficient, 
ks is the Boltzmann constant, T is the temperature and 
At is the timestep. 

It should be emphasized that the purge preserves the 
statistical properties of the unpurged boundary, but will 
not exactly reproduce particle trajectories. Our goal is, 
first, to accurately reproduce the velocity distributions 
at the purge point; and second (and more applicable 
here) , to produce as small a disturbance as possible at the 
purge event. The momentum mirror with strong stochas- 
tic forcing accomplishes both goals. 

Ramping Shock Strength 

The final element of the method is to introduce a qua- 
sistatic change in the shock forcing parameter. The goal 



is to alter the shock strength, while always maintaining 
a direct mechanical coupling between the forcing at the 
piston and the response at the shock interface. The pis- 
ton velocity increases or decreases by a set amount per 
timestep. The point of forcing remains at appoximately 
a constant distance behind the shock front. This process 
continues until the shock strength parameter has reached 
a desired value, at which point the simulation terminates. 

Finally, a second conventional reference run is com- 
pared with the final state of the shock ramp. This al- 
lows the final state obtained from the Continuous Hugo- 
niot Method to be compared directly to a conventional 
shock run to identify any problems and to serve as an 
error check. In our experience, the results of the final 
benchmark runs have been indistinguishable from those 
obtained from the Continuous Hugoniot Method. 

There are two conditions required for the validity of 
our method: first, that the shock must be steady; and 
second, that the state at the purge point must be in 
equilibrium. The second is guaranteed by our thermal- 
ization test before we reduce the system. The first is 
guaranteed, so long as our quasistatic ramp loading is 
slow enough. If the ramp rate is too fast, the system 
is not steady. However, If the ramp rate is quasistatic, 
small changes have time to equilibrate across the entire 
system, and thus the driving and response are mechan- 
ically coupled. If this coupling is not maintained, the 
foundation of the Hugoniot-Rankine equations is eroded, 
and off- Hugoniot states are produced. Our ability to very 
slowly ramp the shock strength parameter while simul- 
taneously assuring that the driving is mechanically cou- 
pled to the shock response is an essential property of this 
method. We accomplish it by driving our system from a 
set distance behind the front. Note that there is only a 
superficial similarity between our method and isentropic 
loading experiements where the driving is applied at an 
ever-increasing distance. If one ramps the shock driving 
velocity in an experimental situation, the result is isen- 
tropic compression rather than a shock response Q, as in 
our technique. 

We can estimate an upper bound for the rate at which 
the shock strength control parameter can be varied. We 
nondimensionalize the velocity by the longitudinal sound 
wave speed behind the shock front Cs , and measure time 
in units of the time needed for sound waves to travel from 
the rear of the system to the front and back again: 
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where Cs is the sound speed in the shocked state behind 
the front. The condition for a quasistatic ramp is then 
given by 



dv p _ d(v p /Cs) 



dt 



< 1 



(7) 
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and thus gives the condition on the velocity ramp rate: 



dv 
~~dt 



c 2 



2L 



dVp Cl 
dt 2L 



(8) 



where the final step substitutes C , the ambient uncom- 
pressed wave speed, for Cs since it is a more easily cal- 
culated quantity and further bounds the value of the ve- 
locity ramp rate. 

Note that the maximum allowed rate that maintains 
mechanical coupling between the forcing and response 
goes to zero for large systems. Thus as L increases one 
not only must deal with increasing numbers of particles, 
but one must also carry out increasingly longer runs. 



600 A in length, consisting of approximately 100,000 par- 
ticles. Continuous Hugoniot Method runs were held to 
200 A in length, consisting of approximately 20,000 par- 
ticles at any one time. The effective system size of these 
treadmilling runs was almost 1.3 //m in length, consist- 
ing of approximately 1,280,000 particles. In all cases the 
shocks were given time to establish a steady state (usu- 
ally 30 to 60 ps). 

In this section, we concentrate our efforts primarily in 
the strong shock (over-driven) regime, v p > 0.75 C Q . This 
single-shock regime presents the fewest complications to 
the co-moving window technique we have developed. Ap- 
plication to double shocks in the elastic-plastic regime is 
discussed in Section V. 



III. APPLICATION TO LENNARD-JONESIUM 

As a first application, we employ the Lennard-Jones 
potential, which has been widely used to study shock 
waves in solids and liquids [1 E El El El El EI 
Although simple in form and application, the Lennard- 
Jones potential exhibits much of the phenomenological 
shock complexity of more "realistic" potentials. We se- 
lected it for two primary reasons. First, there is an ex- 
tensive literature on its shock response, including a pub- 
lished Hugoniot and Hugoniot fit in U s -v p space; second, 
the potential is very fast to compute. 

LJ Simulation Details - We chose to use the cubic- 
spline Lennard-Jones 6-12 potential of Holian EI m 
order to allow easy comparison with published Hugoniot 
results of Germann et al. E3- The spline method uses 
the conventional Lennard-Jones 6-12 form 
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for radii from zero to the inflection point, r sp i, which is 
just beyond the maximum well depth. Between r sp i and 
a point r max , the potential is given by 



= -0-2(1 



r z Y + a 3 (r: 



2 

max 



r 2 ) 3 



(10) 



The conditions for smoothness of the energy and force at 
the inflection point uniquely determine — 0.5424494(7, 
a 3 = 0.09350527ct and r spl = 1.244455cr, r max = 
1.711238cr, where a = 3.405 A. 

The shock was oriented along the (100) direction of 
the face-centered-cubic crystal structure with unit cell 
dimension 5.314 A = 1.561(7. Initial temperature was 
varied with a weak Langevin thermostat from zero to 
10K = 0.083 fcflT/e. Results are found not to depend 
on the initial temperature for shock driving velocities, 
v p above 0.75 C . Systems were 20 x 20 lattice planes 
in cross-section with transverse periodic boundary con- 
ditions. The integration timestep was 0.3 femtoseconds. 

Conventional shock simulation runs (benchmark runs) 
were driven by a warm impactor method and reached 



Piston Velocity Ramp Rate Response - From Sec- 
tion II, the upper bound for quasistatic ramp rate is 



C 2 

v p <^-y = 4.5 x 10 13 m/s 



2L 



= 0.013 m/s per step 



(11) 
(12) 



where the wave speed is given by C Q = ^72e/m [l8| . 
e = 119.8 k B T and m = 6.63 x 10~ 26 kg. L = 200 A is 
the length of the system in these runs. Figure 0] shows 
the results for four ramp rates, on a piston velocity ramp 



from v p 
from v r , 



0.75 C to 1.5 C - The test ramp rates range 
00 (reshock) to 0.077. 
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position (10 meters) 



FIG. 4: Ramp rate determination for Lennard-Jones 

Density profiles for four shock strength ramp rates, ranging 
from reshock to v p — 0.077, are shown. For each rate, the 
density profile evolution is shown solid (v p = 0.75 C ), then 
three later profiles are shown dashed (short dash v p — 0.97 C , 
medium dashed v p — 1.19 C and long dashed v p = 1.41 C ). 

The quasistatic regime is indicated by density profiles 
which are flat and steady. The first three profiles are 
not in the quasistatic regime. The fourth, with non- 
dimensional ramp rate v p = 0.077 (v p = 0.001 m/s per 
step), indicates quasistatic loading for this material and 
system size. 
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The ramp rate of 0.001 m/s per step = 3.3 x 10 12 m/s 2 
(corresponding to v p — 0.077) will be used for the re- 
mainder of the Lennard- Jones work presented. 

A. Principal Hugoniot Results 

The Continuous Hugoniot Method allows a system to 
move directly from shock state to shock state. Therefore, 
the path of our Continuous Hugoniot Method through 
U s -v p (or P-V) space during a single run is the principal 
Hugoniot of hnal shock states in the material. This is true 
to the extent that the system's values within its reduced 
system are indicative of the values very far behind the 
shock (i.e. there is convergence to the final shock state 
within a short distance behind the shock front). Figure 
[S] shows such a path for a run in Lennard- Jones, begin- 
ning at v p = 0.75 C and running continuously through 
Vp = 1.5 C . Initial and terminal reference runs, made 
independently at v p = 0.75 C a and v p — 1.5 C D , respec- 
tively, are also shown. Finally, the Hugoniot fit proposed 
by Germann et al. is plotted over its applicable 
range. 




I i i i i I i i i i I i i i i I i i i i I i I 

0.5 0.75 1 1.25 1.5 

particle velocity (v p /C ) 

FIG. 5: Continuous principal Hugoniot for Lennard- 
Jones - The results of an application of the Continuous Hugo- 
niot Method to a Lennard- Jones system in the range from 
v p — 0.75 Co to v p — 1.5 Co is shown. A plot of the published 
Hugoniot fit of Germann et al. fl^l is plotted in its applicable 
range. The initial and terminal reference runs are also plot- 
ted. Very good agreement is found between the results of our 
method and those of published conventional simulations. 

We see very good agreement of our simulation results 
with both comparisons. In the lower range of v p we can 
compare to the published fit. We see here that our simu- 
lation data overlays the fit very nicely and continues the 
line of the fit beyond the range for which it was origi- 
nally published. At higher values of the piston driving 
velocity we see our simulation data stiffens as it should 
showing a super-linear increase in Us vs u p . In the upper 
range of v p there is no published fit, but we can compare 




position (10" 9 meters) 



FIG. 6: Density comparisons for Lennard-Jones - This 
plot shows density profile snapshots taken from a shock 
strength ramp (black solid) and compares them with density 
profiles from the initial and terminal seed runs (thick dashed) 
and those predicted by the Hugoniot fit of Germann et al. 
(thin dashed). We see near perfect agreement between our 
ramp data and our reference runs. The ramp data is, how- 
ever, consistently below the densities predicted by Germann. 
We attribute this to slow spatial convergence of the density. 

to our reference runs. We can see that the Continuous 
Hugoniot Method data and the terminal reference results 
agree acceptably. 

We would expect that if the Continuous Hugoniot 
Method were following some other path than the Hugo- 
niot, that this would become more apparent at later 
points in the simulation runs (i.e. the error would ac- 
cumulate). If for instance, we were following a thermo- 
dynamic path such as an isentrope the path would di- 
verge and fall below the Hugoniot later in the simulation 
run. This divergence from the Hugoniot states would 
therefore be most obvious late in our simulations. We do 
not see this effect to any significant degree in our data. 
Moreover, we see excellent agreement between the sys- 
tem state of the ramped system and the steady reference 
state at the terminus. 



B. Comparison of Density Profiles 

Figure shows a series of five density profile snapshots 
taken from a single ramp run (shown as solid lines). They 
are equally spaced in the piston velocity from 0.75 C 
to 1.5 C - Also plotted are the density profiles of the 
two reference runs. Using the fit equation and the mass 
conservation jump condition, one can derive the density 
for the 0.75 C state. 

p_ = Us = l + 1.92p p 

Po U s -v p 1 + 0.92v p ^ ' 

For v p = 0.75 C a , the predicted relative density p/po 
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of 1.44 is plotted as a thin dashed line in Figure The 
remaining profiles fall above the applicable range of the 
published fit. 

We see that the density predicted by the published fit is 
something less than 2% above the average density of the 
profiles from the Continuous Hugoniot Method. How- 
ever, in the case of the 0.75 C density profile, we see 
ideal agreement with the seed profile up to 100 A behind 
the front. Beyond 100 A, where the results of our method 
have been purged, the seed density continues to approach 
the fit prediction. We see the same results in a compar- 
ison with the v p = 1.5 C seed run. This observation 
leads us to the conclusion that our method is accurately 
predicting the Hugoniot state density profile within the 
spatial extent of the reduced system. However, densities 
do not necessarily converge within the reduced systems 
to the values predicted far behind the shock front. There- 
fore, comparison between Hugoniot density and the av- 
erage density within our reduced system cannot be made 
directly. If comparison is desired, the density profile of 
the reduced system would need to be extrapolated, or 
the reduced systems would need to be enlarged. 



C. Comparison of Final States 

As previously noted, the final state of the system is a 
particularly important point of comparison because it is 
the state of the maximum integrated error. Direct com- 
parison of every state of the system would be impractical 
due to the computational cost of benchmark runs. We 
can, however, cautiously assume that good agreement of 
the final state indicates good agreement throughout. 




position (10 meters) 



FIG. 7: Temperature comparison - The temperature 
profiles for states obtained from the Continuous Hugoniot 
Method (dashed line) and the terminal seed run (solid line) 
are plotted for shock strength v p = 2000 m/s = 1.5 C in a 
Lennard- Jones system. We see very good agreement in the 
peak temperature and final temperature behind the front. 

We have already seen that the Continuous Hugoniot 
Method produces density profiles which concur with con- 




x position (10 -9 meters) 

FIG. 8: Structural dynamics comparison - Particle slices 
in the region of the shock front for states obtained from (bot- 
tom) the Continuous Hugoniot Method and (top) the termi- 
nal seed run at shock strength, v p = 2000 m/s = 1.5 C in a 
Lennard- Jones system. We see very good qualitative agree- 
ment in the deformation behind the front with pockets of 
order mixed with disorder. We also see the development of 
forward reaching structures ahead of the shock front. 




0.4 0.6 
position (10" 9 meters) 



0.8 1.0 



FIG. 9: Radial distribution function comparison 

RDF for the compressed material extending 10 nm behind 
each front is shown for the Continuous Hugoniot Method 
(dashed line) and the terminal seed run in (solid line) at shock 
strength, v v = 2000 m/s = 1.5 C in a Lennard- Jones system. 



ventional reference runs. Figure shows a comparison 
of the temperature profiles for the final system states. 
Here, the results of the Continuous Hugoniot Method 
are shown dashed and the results of the reference run are 
shown solid, both for a piston velocity v p — 2000 m/s = 
1.5 C Q . The zero position coincides with the shock front 
location. We note, that the temperature peaks are in 
good agreement and the temperature relaxes similarly to 
identical values within the system fluctuations. 

Figure [S] illustrates particle snapshots from the termi- 
nal reference state (top) and the Continuous Hugoniot 
Method (bottom). The slices are along the x-z plane at 
the final piston velocity v p = 2000 m/s = 1.5 C a . The 
particle snapshots clearly illustrate a qualitative agree- 
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ment between the two systems' dynamic structures. We 
see that both systems exhibit a strong disordering transi- 
tion on similar length scales, similar islands of incomplete 
disordering, and similarly sharp density transition. Both 
have also developed forward-reaching features ahead of 
the front. In Figure^ we see the radial distribution func- 
tion for the material behind the front. The rdf results 
for the Continuous Hugoniot Method are shown dashed, 
and the results of the seed run are shown in solid. We 
see from Figure|Hlboth very good qualitative and quanti- 
tative agreement of our method with conventional shock 
methods. 



IV. APPLICATION TO TIN 

This section describes the application of the Continu- 
ous Hugoniot Method to a more realistic material poten- 
tial with the goal of comparing directly with experiments 
on tin. We take advantage of the method's increased effi- 
ciency in producing a dense sampling of Hugoniot shock 
states to study processes at the shock front. 

Tin is a natural choice for laser ablation experiments in 
shock melting because it undergoes shock transitions at 
relatively moderate laser intensities. Experiment places 
its melt on release at ~ 250 kbar and shock melting at 
~ 500 kbar 0]. These pressures are well within the 
range of terawatt class laser systems. Moreover, tin's 
crystal structure makes reflectivity and harmonic diag- 
nostics possible on breakout. 

Modeling Tin with MEAM - We selected the Mod- 
ified Embedded Atom Method (MEAM) as our inter- 
atomic interaction model for tin. This potential has been 
demonstrated to handle shock simulations |20J, and has 
been applied successfully to the study of melt properties 
of tin 21]. 

MEAM was first proposed by Baskes in 1992 as 
an extension of the highly successful Embedded Atom 
Method (EAM) of Daw and Baskes ^ . Both methods, 
inspired by density functional theory, compute energies 
using a semi-empirical combination of two-body inter- 
action and environment-specific electron density embed- 
ding energies. MEAM extends EAM to handle covalent 
bonding by introducing angle-dependent electron densi- 
ties. 

To visualize the embedded atom, recall the chemist's 
notion of the electron density clouds surrounding an 
atom. These lobe-shaped orbital structures overlap as 
atoms approach each other and form a background elec- 
tron density within which each is embedded. With only 
two atoms, this energy of embedding within the other's 
electron density can be absorbed easily into a two-body 
interaction (i.e. bond). As atom density increases, how- 
ever, the electron density becomes a function of many 
atoms and the embedding energy therefore becomes an 
effective environment-dependent interaction. 



The environment-dependent nature of MEAM makes 
it especially good for applications near surfaces, voids, 
defects and interfaces. In these regions, where the lo- 
cal environment is very different from the bulk environ- 
ment, many other potentials are at their weakest, having 
been parameterized with bulk response measurements. 
MEAM adjusts well to these situations, being long-range 
in open structures and short-range in dense closed struc- 
tures, based on screening. 

MEAM has been modified, updated, generalized and 
expanded with subsequent publications |2lL l22l l24j] . Be- 
cause the method varies subtly between publications, we 
specify our implementation in detail in Appendix A. 

Tin Simulation Details - The initial state for our 
runs was a zero-temperature — tin (white tin) in a 
tctrahedral crystal structure with unit cell dimensions 
5.92 A x 5.92 A x 3.232 A l^. The shock was ori- 
ented to run along the <100> direction. All runs were 
10 x 10 lattice planes in cross section. We employed 
periodic boundary conditions in each of the transverse 
directions. The integration timestep was 0.3 femtosec- 
onds. 

Most runs were made using the Continuous Hugoniot 
Method. These runs were held to a reduced system size 
of 165 A in length. Approximately 4,000 particles were 
simulated at any given time. Over the course of the en- 
tire series of runs, the system had an effective length of 
0.26 /im consisting of more than 178,000 particles partic- 
ipating. The velocity ramp runs analyzed in this section 
went from v p = 2000 m/s to 2300 m/s over the equivalent 
of 45 picoseconds (150,000 timesteps). 

Seed runs, using traditional shock simulation methods, 
were driven by a warm impactor. The seed simulations 
were allowed to run for 9 to 12 picoseconds (30,000 to 
40,000 timesteps), until they reached a steady state. Seed 
simulation sizes grew to approximately 580 A in length 
and consisted of approximately 15,000 particles at any 
one time. 

Efficiency of these methods will be addressed in Section 
VI. 

A. Tin Shock Melting and Hugoniot 

Tin undergoes a shock melt and completely thermal- 
izes in simulation on short length scales at moderate 
shock strengths. To verify the properties of the tin melt, 
mean square displacement measurements were made to 
study the temporal dynamics of the melt. This is a very 
good way to distinguish a liquid from a glass or super- 
cooled state. Figure El shows the results of our analysis 
for the same 45 A thick sample as describe above. 

The plot indicates that the material exhibits liquid-like 
diffusional properties. Theory predicts that the mean 
square displacement grows linearly with time in a liquid 

((r(t)-r(0)) 2 ) = 6Dt (14) 
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FIG. 10: Mean squared displacement indicating tin 
melt - The mean squared displacement (msd) is plotted as a 
function of time with a linear fit for a 45 A wide slice of mate- 
rial from behind the front for shock strength, v p = 3000 m/s. 
The linear increase of the msd with time is a classic indicator 
of liquid dynamics. 



where D is the diffusion coefficient. Here D = 
1.84 x 10~ 8 m 2 /s. Experiments by Cahoon found 
a value of D = 2.58 x 10~ 8 m 2 /s for standard pres- 
sure at similar temperatures using an Arrhenius equation 
approach to fit experimental tin diffussion data. High- 
pressure experiments have not been conducted in molten 
tin. 

The velocity distribution functions for the x- 
component of velocity for two moderate shock runs are 
shown in Figure ITT1 Data is shown for v p — 2000 m/s 
(solid line) and v p = 2300 m/s (dashed line). Gaussian 
fits are also shown. These two shock states were chosen, 
as they are the initial and terminal states of a series of 
Continuous Hugoniot Method runs in tin. One should 
notice first that both systems behind the front are well 
described by equilibrium statistics, even in the tails. The 
material has thermalized, a condition necessary for rig- 
orous application of the Continuous Hugoniot Method. 

The fit has the form of a Maxwellian velocity distribu- 
tion 0, 

f(v) = Ae 2 V (15) 

where A is a normalization constant, m — 118.7 amu is 
the particle mass of Sn, and kb is the Boltzmann constant 
and T is the temperature. 

Results in diamonds are for a piston driving velocity 
of v p = 2000 m/s. The fit parameters give values of 
2021 m/s and 3360 K, for v p and T, respectively. Results 
in semi-circles are for a piston driving velocity of v p — 
2300 m/s. The fit parameters there are 2301 m/s and 
4958 K. These are very good fits considering the relatively 
small number of particles available for the calculation. 
Time averaging over 500 timesteps was used to produce 




velocity (m/s) 

FIG. 11: Thermalized velocity distribution functions 

for tin - Velocity distribution functions in v x for two ap- 
proximately 100 A slices of material approximately from 50 
A behind the front. Data is shown for v p = 2000 m/s (solid 
line) and v p = 2300 m/s (dashed line) from a series of runs us- 
ing the Continuous Hugoniot Method in tin. Gaussian fits are 
also shown. The data, plotted in log-linear to accentuate the 
tails, fits nicely, indicating that the material has thermalized 
to a Maxwellian distribution. 



these distributions. 

Recall that we can set an upper limit on the piston 
velocity ramp rate from Equation[S] For tin this gives, 

C 2 

v p < -f = 2.38 x 10 14 m/s 2 (16) 
= 0.0714 m/s per step (17) 

where we have used Co ~ 2800 m/s for the wave speed 
and L = 165 A for the length of the system. This bound 
for the quasistatic piston acceleration is more than five 
times what we derived for Lennard-Jones because the 
wave speeds in tin are higher and the system sizes are 
smaller. Combining this result with the ramp rate testing 
which we conducted in Lennard-Jones, we determine that 
piston velocity ramp rates as high as 0.005 m/s per step 
should be allowed in tin without jeopardizing the qua- 
sistatic condition of the Continuous Hugoniot Method. 

For an added measure of assurance, we use the smaller 
velocity ramp rate of v p = 0.002 m/s per step = 
6.67 x 10 12 m/s 2 in all of the simulations discussed 
in this section. 

Application of the Continuous Hugoniot Method to tin 
gives the path through Us — v p space shown in Figure lT^l 
These results are compiled from a series of runs begin- 
ning at v p = 2000 m/s and running continuously through 
Vp = 2300 m/s. Initial and terminal reference runs, made 
independently, are also plotted. Experimental results 
from the literature are plotted and a linear fit is given 
to the experimental data. The initial shock strength pa- 
rameter of v p = 2000 m/s was selected to assure that the 
results were within the strong shock regime. 
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FIG. 12: Continuous principal Hugoniot for tin - The 

results of an application of the Continuous Hugoniot Method 
to a tin system in the range from v p = 2000 m/s to v p — 
2300 m/s is shown. Experimental data are plotted from sev- 
eral sources and a fit is made to the experimental data. Agree- 
ment is to within 6% across the board, with good agreement 
in slope. 



For experimental comparison, v p = 2000 m/s gives a 
pressure of 730 kbar and v p = 2300 m/s gives a pressure 
of850 kbar. 

We can see that the U g versus v p relationship is ap- 
proximately linear in this range, as in Lennard-Jones. 
The Us response is noisy, but not significantly more noisy 
than traditional shock simulation. Although experimen- 
tal data is not currently available for single-crystal tin, 
we were able to compare to room-temperature polycrys- 
tallinc tin shock experiments. One can see relatively 
good agreement with experimental data from the liter- 
ature j30l l3l| , however our results are consistently above 
experimental values. Agreement with experimental data 
is to within 6% across the board, and shows good agree- 
ment with the slope. 

Given that Ravelo and Baskes found the MEAM po- 
tential to match transition temperatures to within 11%, 
our measurements are likely the best agreement that the 
potential can provide. These results are rather encour- 
aging to further comparison between our computation 
results and experiment. 



B. Continuous Temperature Profiles versus Shock 
Strength 

Our method allows study of the shock front at a fine 
resolution in the shock strength parameter. Where tradi- 
tional simulation and experiment provide discrete Hugo- 
niot points, our method provides a continuum of states. 
This may be particularly useful in the study of melting 
and other continuous transitions. 

We apply the method first to characterize the tempera- 
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FIG. 13: Temperature profile versus shock strength 

for tin - Temperature is represented as a function of both 
position and shock strength. Temperature is represented as 
a third dimension using gray scale. The reference frame is 
moving with the shock front, so that the front is always at 
x = 0. The shock is moving to the right with cold (white) 
pristine material ahead of it. 



ture profiles in tin as a function of shock driving velocity. 
Figure 1131 shows temperature as a function of position 
in horizontal bands, with the shock strength increasing 
up the plot. The shock is moving left to right. The 
white block on the right is the low temperature initial 
state. Temperature profiles were spatially averaged as 
described by Hardy |32| . 

Note first that the temperature behind the front in- 
creases gradually with shock strength. Figure ITU shows 
a slice at a distance 50 A behind the front. We can see 
that the temperature in this region is noisy, but follows 
a clear trend, growing linearly with the shock strength in 
this range. 

More interesting are the temperature profiles within 
the shock front and near the temperature peak. This 
temperature peak is a common feature of the shock in- 
terface. Our data seem to indicate the peak position 
is stable for most of the shock strength ramp. But, at 
approximately 2200 m/s, the position of the peak shifts 
forward abruptly. These results can be seen more di- 
rectly in Figure El We do not have an explanation for 
the behavior. 



C. Melt Length Scale 

The primary comparison point between simulation and 
experiment in this work is the length scale for the melt 
process behind the shock front. Experimental work is 
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FIG. 14: Final Temperature versus shock strength for 

tin - Temperature is plotted as a function of shock strength, 
Up, at a distance of 50 A behind the shock front (hollow). A 
linear fit is provided. Also plotted are temperatures taken at 
constant time (t = .83 ps) behind the shock front (solid) . 
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FIG. 15: Peak temperature position versus shock 
strength for tin - The position, measured from the shock 
front, of the temperature profile peak is plotted as a function 
of shock driving velocity, v p . 
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FIG. 16: Particle melt dynamics comparison - Particle 
snapshots from the region near the shock front where the melt 
is occurring. The length scale for melting can be seen as the 
lattice planes gradually spread to fill the volume uniformly. 
The eye has a difficult time comparing the melt scales. 



cal. Exactly locating a collective response such as a melt 
point is nontrivial, especially with noisy systems. 

Proper quantitative analysis requires calculation of a 
mean square displacement as a function of distance from 
the front. However, no quantitative trend between melt 
length scale and shock strength has been able to be de- 
rived from the analysis, due to excessive noise to signal 
within the system. 




position (10 meters) 



proposed that could measure the timescales between back 
surface arrival of the shock and the melt to within a few 
tens of femtoseconds. This equates to a spatial resolu- 
tion of only a few tens of Angstrom. There is, therefore, 
great interest to establish a relationship between the melt 
length scale and the shock strength, which we have pa- 
rameterized by the piston velocity, v p . 

One can see in Figure Hfll particle profiles for the low 
(v p = 2000 m/s) and high (v p = 2300 m/s) shock 
strength in our ramp. On the right in each is undis- 
turbed crystal and on the left is the shock- induced melt. 
The transition between the two is notoriously difficult 
to characterize, because the nature of the melt is nonlo- 



FIG. 17: Particle collapse - Collapse of all lattice planes 
onto each other indicates the length scale for the diffusion 
of particles to fills space. The Lindemann criterion for the 
onset of 2D melting is met when the mean free displacement 
is greater than a given fraction of the crystal spacing. 

A more visual analysis demonstrates the problem. Fig- 
ure El shows the same data collapsed to a line and av- 
eraged over 1000 timesteps. We can now begin to see 
the length scale to melt developing. If we have in mind 
something like a Lindemann criterion (33 for the onset of 
disorder, we can argue that the melt has occurred when 
a particle from one lattice plane has diffused to a degree 
that it cannot be differentiated from a particle from an 
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adjacent plane. We thereby define the melt point where 
the particles have full coverage of the cell in Figure IT7I 
Thus we can see that the melt occurs on a shortened 
length scale for v p = 2300 m/s by approximately 3 to 
5 A, compared to v p — 2000 m/s. It is clear, now, how 
this signal could be easily swallowed by the noise of the 
system. Further, such a length scale could not be resolved 
exp er iment ally. 

We conclude that the melt front in tin is too abrupt to 
be distinguished from the shock front at shock strengths 
within the strong shock regime. We are challenged to 
move into the weaker elastic-plastic shock regime, where 
the two fronts diverge with time. 

V. TWO-SHOCK STATES AND THE STRONG 
TO ELASTIC-PLASTIC SHOCK TRANSITION 

We begin this section with a brief review of the three 
distinct regimes of shock propagation in simple solids. 
The analysis begins with the shock Hugoniot. For an ex- 
cellent older review of experimental measurements and 
theoretical fits, see Rice, McQueen and Walsh [34j]. They 
provide examples of experimentally determined Hugo- 
niots, EOS models, and the gamut from descriptions of 
experimental techniques to basic theory on shock waves 
in solids. 
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FIG. 18: Shock regime phase space - The three standard 
shock regimes for a solid are identified. The point A is the 
Hugoniot Elastic Limit where the material begins to deform 
plastically due to shear stresses. The point B is the transition 
point on the Hugoniot between the strong shock and elastic- 
plastic shock regimes. 

At high shock strength, the material response is con- 
sidered over-driven. In this case, a single plastic/melt 
wave propagates through the material. The shock front 
shows a steep gradient in all thermodynamic properties. 

As the shock strength is reduced, the material response 
changes and the material supports two waves. There 
is a primary plastic wave front, which is preceded by a 



faster elastic precursor shock. With time, the two waves 
move apart. The division of the shock can be explained 
by looking at the Hugoniot: see Figure EH Points on 
the Hugoniot below the point B are reached by moving 
from the initial condition to point A along a straight line. 
Then a second straight line connects point A to the fi- 
nal state on the Hugoniot curve. One can show from the 
jump conditions in Section I that the slopes of these two 
lines are proportional to the velocities of the two wave 
fronts. In the case of the strong shock (above point B on 
the Hugoniot), a single line connects the initial condition 
and the final Hugoniot state. Thus there is a single front. 

For stresses below the point A in Figure IT81 the solid 
supports completely elastic, reversible shock waves. In 
this regime, transient shock waves which compress and 
release a material leave no permanent distortion. We do 
not address the elastic regime in this work, but assert 
that the Continuous Hugoniot Method could be applied 
to that regime without difficulty. 

Applicability to the two- wave shock - To this point, 
we have investigated only the strong shock regime. Here 
we will extend the technique into the elastic-plastic 
regime and characterize the transition. However, it 
should be made clear that there is some potential dif- 
ficulty in so doing. 

In Section II, we showed that the upper bound for the 
velocity ramp rate used was inversely proportional to the 
system size. This is not a serious hindrance in the strong 
shock regime because we have shown that the system can 
be continually purged, given that the shocked material is 
known to thcrmalizc within a fixed distance behind the 
front. In the two-wave regime we have just discussed, 
however, the elastic precursor can theoretically move ar- 
bitrarily far ahead of the plastic/melt wave. In fact, it 
can be shown that the precursor front must have a higher 
velocity than the plastic front. Thus the system grows 
in time, and long runs must be simulated with large sys- 
tems. 

The assumptions of the method never break down, but 
the computational efficiency advantage of the method 
quickly deteriorates. Larger systems give slower ramp 
rates, which take longer times, which require larger sys- 
tems. We can see from this argument that the method 
must be applied with care to the elastic-plastic regime. 

Results - The results presented here are for a decreasing 
Continuous Hugoniot Method run in tin beginning from 
a shock state just within the strong shock regime. The 
initial state has a shock piston velocity v p = 2000 m/s 
and runs into the elastic-plastic regime to a final state at 
v p = 1250 m/s. 

Figure ^| shows the Hugoniot representation for the 
run. The plot shows that at high shock strength there is 
no precursor ahead of the shock front, but that an elastic 
precursor forms at lower forcing. The transition point 
comes at approximately 1560 m/s in piston velocity. 
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FIG. 19: Tin shock and precursor velocities - Shock 
Hugoniot plot for a decreasing Continuous Hugoniot Method 
run from 2000 m/s to 1250 m/s which shows the emergence 
of the elastic precursor wave. 



We can recast the data by defining an order parame- 
ter as the difference between the shock front and elastic 
precursor velocities. 

$ = Us,precursor ~ Us (18) 

and the nondimensional bifurcation parameter 

g= ^p 1 (19) 

"p, critical 



which seems the most likely conventional scenario to de- 
cribe this noisy system. The bifurcation shows the classic 
linear increase from zero at onset. The bifurcation is con- 
tinuous and imperfect. 

Analysis of fluctuations in approaching onset from 
above shows a characteristic increase in fluctuation mag- 
nitude. Figure shows the standard deviation of the 
order parameter above onset. The inherently noisy na- 
ture of the system makes it possible to fit either a linear 
or sublinear function to the data. 
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FIG. 21: Fluctuation at onset - The fluctuation in the or- 
der parameter increases as the system approaches onset from 
above. 
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FIG. 20: Shock regime bifurcation - Shock regime or- 
der parameter, $ = Us,precursor — Us vs. the shock control 
parameter. The data is consistent with a transcritical bifur- 
cation. 



Figure |23 shows the results of the transformation. 
Solid and dotted lines are overlayed to represent the sta- 
ble and unstable branches of a transcritical bifurcation, 



VI. EFFICIENCY AND SPEED UP 

A rigorous comparison of processor timing between 
conventional methods and our Continuous Hugoniot 
Method is perhaps not possible, because the nature of 
the output data is different. Our method gives a much 
more finely spaced collection of shock states than con- 
ventional methods. 

In the Lennard- Jones system, we saw the computation 
time to compute the two reference runs by conventional 
methods was approximately equal to the computation 
time necessary to compute the entire intermediate Hugo- 
niot, via the Continuous Hugoniot Method. This can be 
a significant speed-up when a dense sampling of states 
is necessary. The impact of the method is also felt in 
memory and disk resources where we were able to simu- 
late the cumulative effect of 1,280,000 particles with the 
resources required to hold only 20,000 at any one time. 

The complexity of the MEAM potential is substan- 
tially greater than that of Lennard-Jones. This complex- 
ity translates into added computation time as processors 
deal with the additional instructions and the multi-tier 
cascade of nested loops. For moderate simulation sizes 
such as 10 x 10 x 40 lattice planes, MEAM has been 
clocked as running between 90 and 200 times slower than 
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a comparable system of Lennard-Jones interactions. The 
factor is significant since it increases shock run compu- 
tational wall clock times from days to months or years. 
The payoff is, of course, that MEAM is a relatively ac- 
curate model for a real material. As such, results can be 
compared to experiment. 

Without the Continuous Hugoniot Method, progress 
in tin might have been prohibitively slow. The data runs 
necessary for a total ramp of Av p = 300 m/s took approx- 
imately 1400 processor-hours on a quad Opteron-64 pro- 
cessor machine using the Continuous Hugoniot Method. 
A conservative estimate for the same work using con- 
ventional methods puts the processing time at approxi- 
mately 20 processor-years on the same equipment. Even 
on a massively parallel supercomputer, the computation 
time is extremely large. 



VII. CONCLUSIONS 

This paper presents and validates the Continuous 
Hugoniot Method by direct comparison of its results with 
published data and with the results of conventional shock 
generation methods. The method was applied to the very 
well studied Lennard-Jonesium test material and to the 
more realistic MEAM potential for tin. 

We found very good agreement in each of our tests 
of the method. We confirmed that the loading path of 
the Continuous Hugoniot Method in Lennard-Jones fol- 
lowed the published Hugoniot fit for both increasing and 
decreasing velocity ramps between v p — 1000 m/s = 
0.75 C and v p = 2000 m/s = 1.5 C . In the range 
where the fit was not applicable, we were able to com- 
pare to conventional simulations methods with equally 
good agreement. By comparing the results of the increas- 
ing and decreasing velocity ramps, we found that the 
system did exhibit some residual rate dependence, but 
was not hysteretic. We found that the method also did 
relatively well in predicting the Hugoniot in the elastic- 
plastic regime. 

Comparison of density profiles produced by our 
method showed very good agreement with conventional 
runs. However, we found that the average density within 
our reduced system is low, on average, by approximately 
2% compared to the final density predicted in the Hugo- 
niot relations. Larger systems would give better agree- 
ment. However, although not converged to the values 
predicted for far behind the shock front, our density pro- 
files do make reliable predictions near the front. 

Comparison of temperature profiles, particle snapshots 
and radial distribution functions at the final state of the 
Continuous Hugoniot Method ramp also demonstrated 
good agreement with conventional shock methods. 

Application of the Continuous Hugoniot Method to the 
study of shock melting in tin clearly demonstrated its 
efficiency and effectiveness. We found that tin both melts 
and comes to equilibrium within just a few nanometers of 
the shock front. Mean square displacements versus time 



showed true liquid diffusion properties in the melt. 

We presented a Hugoniot plot describing the shock re- 
sponse in tin and compared it with experimental data. 
We found that the Hugoniot values for Us versus v p pre- 
dicted by MEAM were approximately 6% above those 
found in experiment. However, differences between the 
experimental parameters and the simulation runs could 
account for such a discrepancy. We plotted the temper- 
ature profiles of tin as a function of the shock strength 
parameter, v p and noted an anomalous transition at high 
driving velocity in which the temperature peak migrates 
abruptly forward in the shock profile. 

The Continuous Hugoniot Method was applied to 
shocks in the two-shock elastic-plastic regime. The tran- 
sition from the strong shock regime was shown to be 
continuous and consistent with a transcritical bifurcation 
with order parameter equal to the relative velocity of the 
shock and precursor waves. The fluctuations in the or- 
der parameter grew as the system approached onset, as 
expected. 

Finally, we were able to make these measurements with 
greatly reduced computational expenditure over conven- 
tional methods. These savings proved critical when the 
method was applied to the more realistic and computa- 
tionally costly MEAM potential. 
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APPENDIX A: MEAM POTENTIAL 

The total energy given by MEAM is 



E = H kfei + jEw ( Al ) 

* \ j¥» / 

where the sums are over particle indices; F is the embed- 
ding energy as a function of pi, the background electron 
density at the site of the i th particle; and is the two- 
body interaction between particles % and j as a function 
of Rij, their separation distance. 

The embedding function, F, has the form 

F{pi) = AE c pi\nPi (A2) 

where A is a free parameter and E c is the cohesive energy. 
The background electron density for tin is calculated at 
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each occupied site by 



where p^ is the partial electron density associated with 
the spherically symmetric s orbital contributions from 
surrounding atoms. 

r is a weighted sum of the non-spherically symmetric 
partial electron densities associated with the p, f and g 
orbitals. T is positive definite for tin. 

h=l VP 7 

ft 1 ', t( 2 \ and are parameters indicating the relative 
importance of each orbital, and the higher partial elec- 
tron densities. 

The second term from the MEAM total energy, Equa- 
tion ^d] is the two-body interaction term 

ct>(R) = l{E-(R)-F(p (R))} (A5) 

where Z is the number of nearest neighbors in the refer- 
ence structure and F(p°(R)) is the embedding energy of 



the reference structure background electron density p°. 
E U (R) is the energy per atom in the reference structure 
as a function of the nearest-neighbor distance R. 

In lieu of a potential cutoff distance, MEAM imple- 
ments a many-body screening. Thus the effective cut- 
off distance depends on the local environment. Dense 
structures have short cutoffs; and open sparse structures, 
where there is little screening, can have long-range inter- 
actions. 

The screening function < < 1 multiplies the elec- 
tron densities and the pair potential. The total screening 
function is the product of screening terms for all particles 
j which reside between the i and j particles. The degree 
of screening is determined by 

Cik = J\ Sl 3 k ( A6 ) 

where Sijk is the screening effect of j between the i th and 
j th particles. We use the simple algebraic form [2{| for 
the screening term. 

We take our MEAM parameters from Ravelo and 
Baskes |21l without modification. 
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